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ABSTRACT 

An acoustic radiation analysis has been developed to predict the 
far-fleid characteristics of fan noise radiated from an 
acoustically lined unflanged duct. This analysis is comprised of 
three modular digital computer programs which together 
provide a capability of accounting for the impedance mismatch 
at the duct exit plane. This report discusses the Duct 
Termination Impedance Program whose relationship with the 
other two modular reports of the analyses is illustrated on the 
following page. 

Admissible duct configurations include circular or annular, with 
or without an extended centerbody. This variation in duct 
configurations provides a capabiPty of modeling inlet and fan 
duct noise radiation. 
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DEVELOPMENT OF ACOUSTIC RADIATION ANALYSIS 
OF TURBOFAN NOISE FROM LINED, UNFLANGED DUCTS 
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1.0 


INTRODUCTION 


1.1 Background 

This report discusses the Duct Termination Impedance Program that forms the 
intermediate link for the package of modular computer programs which together form a 
first generation capability for the analytical prediction of far-field noise being radiated 
from an unflanged inlet or exhaust duct. 

The development of this program had its origin in the need for a method to correct the 
far-field directivity patterns for the discontinuity in the duct impedance that exists at 
the duct termination. It is known that when a propagating wave traverses the duct exit 
plane that part of the energy Is reflected back up the duct while the remainder is 
transmitted through the interface. In addition, the radial variation of the termination 
impedance generates higher radial duct modes. Thus, both reflected and transmitted 
waves contain modes other than the modes comprising the incident wave. Obviously, it is 
necessary to understand the concept of modal coupling to determine the amount of 
energy radiated to the far-field and its modal content 

The subject program couples the infinite space exterior to the duct, to the interior of the 
duct by accounting for the impedance mismatch at the duct exit plane. This digital 
orogram utilizes the output from the Acoustic Source Distribution Program, Reference 
1, to generate the pressure distribution and impedance which act on the duct exit plane. 
This information is then used by the Directivity Index Program, Reference 2, to 
determine the far-field radiation patterns for a semi-infinite length duct which are 
corrected for the impedance mismatch at the duct exit plane. 


1.2 Technical Approach 

This technical development is concerned with the effect of the impedance mismatch at 
the duct exit plane on the radiation of fan or compressor noise from an inlet or exhaust 
duct. Admissible duct configurations, shown in Figure 1, include unflanged circular or 
annular ducts, with or without extended centerbodies. 

The approach presented here is an extension of the case of the radiation of sound from a 
duct with an infinite baffle on its end to the unflanged case which accounts for the 
diffraction of sound around the outside lip of the duct. Phased contributions of the 
outer and extended inner duct surfaces and the duct exit plane must be determined and 
then combined to determine the total pressure acting on the open face. 

The duct termination impedance characteristics of the duct configuration are determined 
from the spatial distributions of acoustic sources over the duct walls and exit plane. The 
source distribution which satisfies the imposed boundary conditions is determined by 
the Acoustic Source Distribution Program, Reference 1. These boundary conditions 
require that the acoustic perturbation velocity normal to the outer duct walls must 
vanish. In addition, the normal velocity distribution over the duct exit plane must equal 
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the velocity distribution associated with either a single mode or a combination of modes 
of discrete frequency. 

The acoustic source strengths are numerically integrated over the duct walls and exit 
plane to yield the pressure distribution on the duct face. From the calculated pressure 
distribution and the known velocity distribution, two impedances are calculated. The 
first is a localized or point impedance which is defined as the ratio of the pressure on a 
given segment on the duct face to the velocity on that segment. This distributed or local 
impedance reflects the modal coupling effects of the impedance mismatch at the open 
end of the duct. The second impedance is a modal impedance which is constant across 
the duct face and provides a measure of the radiation efficiency of the given velocity 
distribution. 

Several assumptions have been made in this development. These concern the transport 
characteristics of the model and the decomposition of the acoustic sound energy into the 
acoustic modes present in the duct. The acoustic sound is assumed to be radiated 
through an ideal fluid (nonviscous) which has no mean flow. This simplification will not 
significantly affect the results for noise propagation through a locally subsonic flow field 
of an inlet. However, radiation through an exhaust or jet flow will not be represented 
properly by this model. 

The velocity distribution Is assumed to have a Cos (m0) type of angular dependence. 
This simplification is justified by the angular characteristic functions determined from 
the solution of the governing wave equation for sound propagation through annular 
ducts. This assumption does not affect the ability to couple different modes in the 
far-field program, Reference 2, but does restrict this program to the extent that for a 
velocity distribution whose angular dependence is specified by a combination of 
different m angular modes, a separate analysis must be made for each set of angular 
modes. 
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2.0 


MATHEMATICAL DEVELOPMENT 


The acoustic pressure on the open face of an unflanged duct is formulated in terms of 
the acoustic velocity potential. This formulation provides a mathematical model with a 
straightforward method of satisfying the velocity boundary conditions on the duct walls 
and exit plane. The phase of the pressure wave transmitted through the duct exit plane is 
accounted for in the phase relationship of the potential function. 


Since the velocity potential is governed by a linear set of equations, superposition may 
be used to determine the field potential 0 in the region exterior to the duct. The 
potential for the mathematical model shown in Figure 2, is determined by the 
superposition of three potentials 0 q, ^ , <f> 2 , which represent the contribution of the 
duct face, outer duct wall and extended centerbody duct wall, respectively. This 
technique is expressed mathematically by 

0 (r t =? + (1} 


The component velocity potentials are determined by a spatial distribution of monopole 
acoustic sources on each reflecting surface. For example, the velocity potential due to a 
distribution of harmonically vibrating sources on the face of the duct is 
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The form of this integral equation is derived from the following assumptions; 

• Acoustic velocity at the point P may be determined from a monopole 

acoustic source kR 0 )^. w h\ c h represents an outward 

propagating wave and is a solution of the wave equation in spherical 
coordinates. 

• The angular dependence of the potential may be described in terms of 
the characteristic function Cos (mf?) for a given angular mode number m. 

• The radial distribution of sources is given by the spatial distribution 
function 3> 0 (r} which is determined by the Acoustic Source Distribution 
Program, Reference 1. 


Similarly, the velocity potential terms for the outer and extended inner duct wall 
surfaces, respectively, are 
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Acoustic pressure is related to the velocity potential by 

d<P 


P = P 


dt 


The velocity potential of Equation (1) yields the following expression for pressure 
varying harmonically. 


p = — LcoP | <4+ 0 2 1 


(5) 


To determine the pressure at the point r ,6 on the open face of the duct, the separation 
distances in Equations (2) through {4} are 
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where j = 0, 1, 2 


^o r o 
A-j = b 

A 2 = a 


A useful quantity in the study of sound radiation from ducts is the modal impedance in 
the axial direction. The modal impedance at the duct exit plane is defined in terms of 
the acoustic pressure and velocity, in any mode, which are present on the duct exit 
plane. 

z=-H- 

v 

» 

The modal impedance which is expressed mathematically as a complex quantity, is 
written in terms of its real and imaginary parts which are normalized by the specific 
impedance of ai r pc, by 

Z= Pc(t-lX) j 7 j 

The variable r is the radiation resistance which is a measure of acoustic energy radiated 
through the exit plane and the variable X is the radiation reactance which is a measure of 
the acoustic energy reflected back up the duct. 

Acoustic pressure and velocity are related to the impedance through the sound intensity 
radiated through the duct exit plane. In a stationary medium, the sound intensity is 
defined to be the time average of pressure times velocity. Mathematically, the intensity is 
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related to the real part of the product of the complex conjugate of the pressure and 
velocity by 
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Sound intensity, normal to an area, is the rate at which the sound energy crosses a unit 
of that area. Thus, the total sound power radiated through the duct exit is 
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Sound intensity is related to the impedance by two expressions; one which involves only 
the radiation resistance 
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and the other which includes the radiation reactance. 
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These expressions are substituted into Equation {9) and the resulting expressions for 
radiated energy are equated to yield the following relations for the impedance: 
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METHOD OF SOLUTION 


The constraint on this program to be compatible with the Acoustic Source Distribution 
Program of Reference 1, required that the "Box Method" technique be used to solve 
Equations (2) through {4}. This mathematical technique assumes that the variation of 
the source distribution over a small portion or box of the reflecting surface is constant, 
an assumption which allows the integrations to be evaluated by a summation. To 
facilitate the numerical solution of these equations, the duct exterior face is divided into 
a sequence of M annular rings, and the duct walls into N and N2 cylindrical segments as 
shown in Figure 3. 

Application of the box method to 0 q, the duct face contribution to exit plane pressure, 
enables the integral equation to be written as 

= cu n 

M 

% 

1=1 

For the duct wall contributions, tne equations become 



Equations (12) through (14 ) are integrated numerically. The numerical technique used 
is an even order Gaussian quadrature scheme. The integration for each annular segment is 
performed by an iterative method. An initial even order quadrature is chosen (usually 
2nd or 4th order} and the integral evaluated. The order is increased by 2 and the 
integration repeated. 

The two values of the integral are compared. If they agree to within a specified error, the 
next segment is integrated. If agreement is not obtained, the order is again increased and 
the procedure repeated up to a specified maximum order. At this point, the last 
calculated value of the integral is used and the integration procedure passes to the next 
segment. 
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The quadrature method used is the Gauss formula for arbitrary intervals. 
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This relation states that the integral may be approximated by a summation of NG (the 
order of the quadrature) evaluations of the function at a prescribed abscissa Xj which are 
multiplied by the associated weight Wj. The abscissae and weights of the Gaussian 
Quadrature may be found in Reference 3. The referenced abscissae have the normalized 
range -1 < y < 1 and are related to the arbitrary interval abscissae by 

(b + a) (b -a) 

x j= — 

The numerical equations for the evaluation of the contributions of the duct face, outer 
duct wall and extended centerbody wall to the pressure, denoted by pg, p.| P 2 , 
respectively, are 
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where the quadrature points are 
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operations of Equation (5) for each annular ring. 

Pi = P'(Pq + P-j + P2h 


8 


In the evaluation of p 0 , the integral cos (^R 0 )/R 0 in Equation (16) is nearly singular 
when the integration is performed over the box containing the control point. The effect 
of this singularity on the numerical quadrature is a degradation in its convergence. If the 
error resulting from a low order quadrature is not acceptable, either a higher order 
quadrature or a closed form evaluation of Equation (16) must be used. 


In view of the extremely poor convergence of the quadrature in Equation (16), a closed 
form approximation of the integral was used. A technique was developed that consisted 
of removing the singularity by a change of variables and using the series expansion for 
the exponential function. The excellent convergence properties of this approximation 
required only a small number of terms. This integral equation is approximated by 
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In is tlie natural logarithm 
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^ = 1 - (WA )Vl2 + (WA)V320 - (WA) u /l6l2S 
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fj. 0 - 1/6 - (WA) 2 (sec 2 a +3/2)/24o+(WA)^(see\x_ 5_sec 2 a +5/8)/504o 
* 3 12 

^ = 1 - (WB) 2 /12 + (WB)V320 - (WB) 6 /l6l28 
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The local impedance for each annular ring Is calculated from the acoustic velocity and 
pressure determined at the duct termination. 


z = A 
1 v. 


The radiation resistance and the radiation reactance of the modal impedance are 
calculated from the numerical form of Equations (10) and (11). 
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4.0 


RESULTS 


4.1 Application 

The numerical techniques developed in Section 3 are applied to the mathematical 
representation of the radiation problem developed in the Acoustic Source Distribution 
Program, Reference 1. Duct segmentation parameters and corresponding source 
distributions are Input to this program to determine the termination impedance 
characteristics of the duct configuration. Far-fieid directivity patterns corrected for 
termination impedance effects are determined for each mode and then combined to 
yield the total noise signature. In the event both an inlet and fan duct are to be modeled, 
their respective far-field directivity patterns can be calculated and then combined to 
yield the total noise signature of the nacelle. 

An example duct radiation problem is presented to illustrate the output of this program. 
A source distribution which was generated by Reference 1 for a plane wave velocity 
distribution on the face of a circular duct (example problem 1) was input to this 
program. The resulting pressure distribution for the duct exit plane is shown in Figure 4. 
Value: of pressure which are normalized by the characteristic impedance of air, pc are 
plotted against the radial position of each control box. Behavior of the pressure near the 
outer duct wall is due to the singular behavior of the source distributions, which is 
discussed in Reference 1. Since the velocity distribution for a plane wave is uniform 
across the duct, the resulting local impedance is proportional to the pressure. Thus, the 
local impedance is shown to vary across the duct cross section. The value of modal 
impedance that resulted from this analysis is given by its components: 

T = ,520 

X=1.93 


4.2 Limitations 

There are several limitations to this analytical program. First, the effect of mean flow 
within the duct on the termination impedance has been neglected. Although application 
of this analysis to a nacelle radiation problem does not yield an exact model, it is 
anticipated that this simplification will not significantly affect results for noise 
propagation through a locally subsonic flow field of an inlet. However, radiation through 
an exhaust or jet flow will not be represented properly by the present analysis. 

The second limitation to the program is restriction of the specified velocity distribution 
on the face of the duct to a Cos (m<?) type of angular dependence. This assumption does 
not affect the capability of the program to combine different m angular modes in the 
far-field by superposition of the source distributions. It does, however, restrict this 
program to the extent that for a velocity distribution whose angular dependence is 
specified by a combination of different m angular modes a separate analysis must be 
made for each set of angular modes. 
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5.0 


CONCLUSIONS AND RECOMMENDATIONS 


An acoustic radiation analysis has been developed that predicts the angular distribution 
of acoustic energy in the far-field that is radiated from an unflanged duct. In spite of its 
many advantages, this acoustic radiation analysts is limited in its range of application. 
For this reason, this analysis has been relegated to the role of a first generation analysis 
in a commitment to develop an analytical package which will predict radiation 
characteristics of a turbofan jet engine. As discussed in Reference 1, the present version 
of the analysis reaches optimum efficiency when it is applied to low frequency noise 
radiating from an infinite duct This and the other imitation due to a Cos(m0) type of 
angular dependance of the velocity distribution are an inherent property of the "box 
method" numerical technique. 

A second generation version of the program could avoid these limitations by utilizing a 
collocation technique of evaluating the integrals of Equations (2) through (4). The 
primary feature of this technique is that a selected analytical expression is used to 
represent the source distribution. The time saving feature of the collocation technique is 
that it requires fewer evaluations of the integrals. The collocation procedure would 
enable the analysis to handle any combination of angular modes, a result which would 
remove the angular modal coupling limitation and eliminate the need for a large number 
of angular boxes. 

in addition to the collocation procedure, a second generation version of the analysis 
should contain a method of accounting for a mean flow within the duct. This method 
should account for the effects of velocity gradients on the distortion of the radiation 
field by the jet. 


|BKCEDJNa PAQiJ 
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NOMENCLATURE 


| j 


■ I 

Computer Program Variables ■ ! 


A 

ALIMP 

B 

C(l) 

CHI 

C2(I) 

CHK 

D{IJ 

INK 

K 

L 

MW 

M 

N 

N2 

PO) 

PH 1(1} 
R 

RO,RI 

RHO 

S 

TAU 

TT 

Tl 

VU) 

w 

WDS 

XANG 


Centerbody Radius 
Local Impedance 
Outer Duct Radius 

Dimensionless Size of ith Box on Duct Outer Wail 
Radiation Reactance 

Dimensionless Size of Ith Box on Extended Centerbody Wall 
Convergency Criterion 

Dimensionless Size of Ith Box on Duct Exit Face 
Configuration Parameter 

No. of Polar Angles at which Far Field Data is Required 
No. of Angular Boxes 
Angular Mode No. 

No. of Angular Boxes 

No. of Cylindrical Boxes on Duct Wall 

No. of Cylindrical Boxes on Extended Centerbody Wall 

Pressure Distribution of Ith Box 

Source Strength Distribution of Ith Box 

Radial Coordinate 

Distance between Point on Duct and Field Point 

Density of Air 

Sound Speed in Air 

Radiation Resistance 

Angular Coordinate 

Angular Coordinate Integration Variable 
Velocity Distribution of Ith Box on Exit Face 
Frequency Radians/Second 
Wave Number W/S 

Nondimensional Angular Size for Control Box 
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The basic operation of the Duct Termination Impedance Computer Program is described and the 
input/output formats are listed. A flaw chart of the program is presented in Appendix 2 and listings 
of the main program and its subroutines are presented in Appendices 3 through 7. An example 
problem which illustrates the input preparation and output format of the program is presented in 
Appendix S. 


INPUT 

The input data for the program consists of 23 variables whose symbols and card formats are 

described below: 

CARD OR 

CARD SET FORMAT DATA 

1 3612 K Number of angular points, at which the directivity 

index will be calculated. 

M Number of annular rings used to represent the face of 

the duct. 

N Number of cylindrical rings used to represent the 

outer duct wall. 

N2 Number of cylindrical rings used to represent the 

extended inner duct wall (if present) 

INK Indicator for extending the inner duct wall past the 

duct face. 0 does not extend the wall, 1 extends it to 
infinity. 

MW Angular mode number. 

NN Maximum number of passes to be made through the 

quadrature iteration loop. NN is related to the highest 
quadrature order NG by NN = NG - 2 . 

2 

L Number of angular segments. 

LP Number of integrations in integration routine. 


2 6E12.6 A Inner duct radius 

B Outer duct radius 

W Circular frequency, radians/second 
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S Speed of sound in units consistent with those of A 

and B 

CHK Convergence criteria for numerical integration 

XANG Nondimensional angular size for control box Radians 

= 2 jt(XANG}. 

RHO Air density in units consistent with A, B, and S. 

3 6E12.6 C{l), Outer duct wall box lengths. These lengths are 

I = 1, N normalized by the outer duct radius, numbered with 

C(1) closest to the duct exit plane and I increasing in 
the negative z direction. 

4 6E12.6 DO), Annular ring box widths. These widths are 

l = 1, M normalized by the outer duct radius, numbered D(1) 

at the duct centerline (or adjacent to the centerbody 
wall duct exit plane intersection/and I increasing to 
D(M) at the lip of the duct. 

5* 6E12.6 C2(l} Extended centerbody wall box lengths. These lengths 

l = 1, N2 are normalized by the outer duct radius, numbered 

with C2(l} closest to the duct exit plane and I 
increasing in the positive z direction. 

6 6E12.6 V(l), Velocity distribution on the duct exit plane. This 

I — 1, M distribution is numbered with V(1) at the duct 

centerline (or adj'acent to the centerbody wall duct 
exit plane intersection) and I increasing to V(M) at 
the lip of the duct. 

7 6E12.6 PH 10 Source distribution on duct face. Method of 

[1) = 1,M numbering this source distribution ts analogous to 

D{1). 

8 6E12.6 PH110), Source distribution on outer duct wall. Method of 

I = 1, N numbering this source distribution is analogous to 

C(l). 

9* 6E12.6 PH 12(1 ) Source distribution on extended centerbody wall. 

I = 1, N2 Method of numbering this source distribution is 

analogous to C2(l). 

*Card sets 5 and 9 are needed only if INK = 1. 
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PROGRAM FLOW 


Program operation consists of two separate parts. The first part is a loop which calculates the local 
pressure distribution and the resulting local impedance. Equations (IS) through (IS) are solved to 
yield the exit plane and duct wall contributions p Q , p^, P 2 to the pressure distribution on the duct 
exit plane. The second part calculates the value of the modal impedance. Equations (19) and (20) 
are solved to yield the resistive and reactive components of the modal impedance. 

Four subroutines are required to calculate the pressure distribution on the duct face. PRESO, 
PRES1, ar .;2 PRES2 calculate the contribution of exit plane and duct wall contributions to the 
pressure distribution. GAUSS is the subroutine that provides the weights and abscissae for the 
Gaussian quadrature described in Equation (15). Listings of PRESO, PRES1, PRES2, and GAUSS 
are presented in Appendices 4 through 7, respectively. 

OUTPUT 

The output format consists of three groups of data: 

1. Input Data — All pertinent input data are printed out. 

2. Distribution of pressure and impedance — Values of local impedance, pressure and velocity 
for each annular ring and its corresponding radial position are printed out in row format. 
For each angular integration that fails to meet the convergence criteria, values of angular 
position, radial position, differences of the real and imaginary parts of the integral are 
printed out. 

3. Modal Impedance — The resistive and reactive parts of the modal impedance are printed. 
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APPENDIX 3 


MAIN PROGRAM LISTING 


C THIS PROGRAM CALCULATES THE ACOUSTIC IMPEDANCE OF 

C THE DUCT EXIT PLANE 

C 7DM915 20 - 20 

INTEGER *4 DRECTP/02/ 

COMPLEX *8 CMPLX 

1 , PO ,P 1 , P2 ,P 

2 ,P H lOt 40) t PHI1 (40) ,PHT 2 140) ,V(40),ALIMP 
DIMENSION 

1 PN( 40 ) , PD ( 2 » 1 ) * VD (2 ,1 ) , VMC4C ) , PM ( 40 ) 

2 ,D(40) ,C (40) ,P (80) ,C?(40) f E(40) 

EQUIVALENCE (PD(1,1 ),P{ I) ) ,IVDUfl),V(I)) 

COMMON /CAT A/L,M,N,NN,LP ,MW, A,B,W,S,CHK,D,t ,E ,PI 
COMMON /PHI /PHI 0» PHI 1, PH l 2 , PO, PI ?P2 

COMMON /C2N2/ C2 f N2 

COMMON/ANG/ FINK, BD2 , WDS, FMW, PI2, BW, AW ,8B,AA 

100 FORMAT 0612) 

220 FDRMAT (6E12..5I 

10C0 FORMAT Cl TUBE END IMPEDANCE PROGRAM FOR UNFLANGED ANNULAR DUCT 
1 *' // 25 X, 'INPUT DATA') 

402 FORMAT ('0 THETA BOX WIDTH ON DUCT FACE'/) 

403 FORMAT Cl VELOCITY ON DUCT FACE'/) 

407 FDRMAT Cl*) 

4Q8 FORMAT CO A = *, 1PE13.6, « B = *, E13.6, ' W =', E13*6, * S = ', 
1 E 1 3* 6/ * OC HK- 1 , E13.6, ' RHO=' , E13.6) 

1100 FDRMAT CO ANNULAR RING BOX WIDTHS* /) 

1200 FDRMAT CO OUTER DUCT WALL BOX LENGTHS' /) 

1300 FORMAT CO INNER DUCT WALL BOX LENGTHS' /) 

1400 FDRMAT CO SOURCE DISTRIBUTION ON DUCT FACE* ' /) 

1500 FORMAT ('0 SOURCE DISTRIBUTION ON OUTER DUCT WALL' /) 

1600 FORMAT CO SOURCF DISTRIBUTION ON INNER DUCT WALL* /) 

C * * * * # * * ****** *** 

READ (5,100) K, H, N, N2 , INK, MW, NN, L,LP 

READ ( 5,220) A, B, W, S, CHK , XANG , R.HO 

READ (5,220) (C(I), 1=1, N) 

READ (5,220) (0(1), 1=1, M) 

IF (INK „GT* 0) READ (5,220) (C2(I )"I=1,N2) 

READ (5,220) (VIII, 1=1, M) 

READ (5,220) ( PH 1 0 (I) , I =1 , M) 

RE AO (5,220 ( PH 1 1 (I) , I = 1 , N) 

IF (INK .GT. C) READ (5 ,220) (PHI 2 ( I I ,1 = 1, N2) 

E{ 1) = XANG 

DT2= l l.OEO - XANG) / (L-i) 

DO 50 1=2, L 
50 Ft I I = DT? 

10 5 FORMAT (1PE16.6, 3E18.6) 

104 FORMAT CC L =', T4, * M =*, 14, » N =', 14, • NN =• , 14, 

1 ' LP 14, ' MW =*, T 5 , * INK =», 15/) 

WRITE (6,1000) 

WR ITE( 6, 408 ) A,B,W,S,CHK, RHO 
WR I TE ( 6 , 104 ) L,M,N,NN,LP,MW,INK 
WRITE(6,105) (DM), 1=1, M) 
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WRITE {6, 1200) 

WRITE{6,105) (CU)*I«lfN) 

IF (INK .GT. 0) WRITE (6, 13001 

IF (INK * GT* 0) WRtTE (6, 105) (C2U)t I = lvN) 

WRITE <6, 402) 

WRITE(6,105) (Ell), 1=1, L) 

WRITE (fit 403) 

WR.ITECfif 1C51 (V(I),I=1,H) 

WRITE (fit 1400 i 

WR ITE ( 6 , 105 ) IPHIOC II*I«1»MJ 

WRITE (fit 1500) 

WRITE! 6,1051 ( PHIll I) , I -1,N> 

IF l INK .GT. 0) WRITE. (6, 1600) 

IF (INK .GT. 0) WRITE (6,105) ( PHI2 (I) , 1 = 1 ,N2 ) 

C** INITIALIZATION * * * * 

PI*3. 141592653589793 
FINK = FLQATl INK) 

B02 = B / 2. DC 
WDS = W / S 
FMW = FLOAT(HW) 

P 1 2 * PI * 2. DO 
BW = B * W 
AW = A * W 
AA = A * A 
BB = B * B 

THI5 LOOP CALCULATES THE RADIAL PRESSURE DISTRIBUTION 

DD 200 JC = 1, M 

P(JC)= CMPLX( 0. CEO, 0. OE 0 ) 

IF t JC-1 )210,210 T 211 

210 R=DUC)*B02 LA* FINK 
GO TO 212 

211 JE=JC- 1 

R=R+(Dl J E ) +D3 JC ) )*BD2 

212 CONTINUE 

IF CR.GT.AI GO TO 170 

IF (INK.EQ.l .OR. VD ( 1 , JC ) * EQ. 0. DO ) GO TO 203 

151 FORMAT! 1HC,5X, 13, 6X,D25. 14) 

170 CONTINUE 

CALL PRESO(R) 

CALL PRES HR ) 

IFl INK-1) 400, 5G0t 500 
500 CALL PRES2IR ) 

P( JC)= (P-0 4 PI *- P 2) * RHO 
GO TO 201 

400 P( JC) - IPO + PI) * RHO ' 

201 CONTINUE 

ALTMP = P( JC) / V( JC) 

GO TO 202 

203 P( JC>= CMPLXt 0. CEO, O.OEO) 



202 CONTINUE 

WRITE (6, 160) JC, Rt ALIMP, P(JC), VIJC) 

160 FORMAT (* ANNULAR RING NO*, 13, 1 RADIAL POSITION =',1PE14.6/ 

1 25X, ’LOCAL IMPEDANCE = ' » 2E14*6/ 

2 32X , ’PRESSURE = ', 2E14. 6/ 

3 3?X, ’VELOCITY = ', 2E14. 6/) 

20C CONTINUE 

C 

C THIS LOOP INTEGRATES THE VELOCITY AND PRESSURE DISTRIBUTIONS 

C OVER THE DUCT FACE 

C 

DO 300 1=1, M 

VMM)=VDt 1, 1)**2+VDI2,n**2 
PN( l)=PD< 1. I )**2+P0 (2,1 )**2 
300 PM t T ) — P D C 1, I )# VDt I, I ) + PD12,I) * VD(2rD 
RINT=0. C 
P INT=0 • C 
VINT=O.C 
DO 330 K- 1 , M 
IF ( K- 1 ) 331,331, 33? 

331 RR=D( K )*BD2 + FINK * A 

GO TO 333 

332 KL=K— 1 

RR=RR+ ( D( KL ) +0 { K ) J+BD2 

333 CONTINUE 
RA=RR-0{K)*BD2 
RR=RR+D(K)*8D2 

DEL= ( RB*RB-RA*R A) /2.Q 
R INT=RTNT+PNtK)*DEL 
PINT=P INT+PM( K )*OFL 
330 VINT=VINT+VM(K )*DEL 
C 

C CALCULATE THE REAL AND IMAGINARY PARTS OF THF 

C TUBE END IMPEDANCE 

C 

T4U=P IN T / 1 VINT* S+RHO) 

C H I = ( ( VINT*RINT-PINT*PINT)**0,5) / ( S*VI NT*RHO ) 

WRITE16,305) TAU 

306 FORMATt 1H , 10X, ’RESISTA NCE RATIO, TAU= **£15.8) 

WRITE( 6,3050) CHI 

3050 FORMAT ( 1H , 10X , ’RFACTANCE RATIO, CHI = »,E15.8) 

5000 CONTINUE 
STOP 
END 
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SUBROUTINE PRESCIRJ 

7.0M915 21 - 21 

THIS SUBROUTINE CALCULATES THE CONTRIBUTION OF PHTO 
TO THE PRESSURE ON THE DUCT EXIT PLANE 
REAL *4 F5D8/.625D0/, F 5D12/* 41 6 66 666666666 66D0 / 

1 , F 16D5 /3.2D0/ ,F6D5/l*2DO/ 

COMPLEX *8 CMPLX 

1 ,PHI0I40) t PHI1I40) r PHI2 140) »P0,P1,P2 t PART 

DIMFNS ION X(20) ,WT( 20) , D ( 40 J ,C f 401 T E 140 ) T PB £40 ) , PT( 40 ) 
COMMQN/PHI/PHIQ»P'HIl,PHI2tPO»Pl »P2 
CnMMON/DATA/L,M,N,NN,LP ,MW,A,B,W,S,CHK,D,C ,E,PI 
COMMON/ANG/ FINK, RD2, WDS, FMW, PI2, BW, AW f BB, AA 
PO^ CMPLX(C.OEC,Q.OEO} 

CT = 1 # 0 
ST=0 . 0 

LOOP FOR THE ANGULAR VARIABLE 

DO 220 KC=1»L 
IFCKC-1 1221,221,222 

221 TT = 0 . 0 

GO TO 222 

222 XCL =KC- 1 

TT = TT + ( EtKCL) *E (KC ) )*PI 

223 DT=E t KC )*2.0*PI 

LOOP FOR THE RADIAL VARIABLE 

□0 230 K= 1 »M 
I F ( K— 1 ) 221,231,232 

231 RR=D(K)*BD2 +A* FINK 
GO TO 232 

232 KL=K-1 

RR=RR+( D(KL)+D(K) )*BD2 

233 CONTINUE 

CLOSED FORM INTEGRATION OVER THE PRESSURE SINGULARITY 

BZ=D{K 1*BD2 

A7=R*DT/2.C 

BZDT - 8Z * DT / 2.00 

I F I K C— 1)234,234,235 

234 RQ=R+DI KJ*B/5.0 
RP=R-D( K )*B/5.0 

[Ft {RR .GE.RQ ) .OR. (RR* Lt «RP! 1 GO TO 235 

TN=BZ/AZ 

AL P= ATAN(TN) 

SEC=1. 0 / COS{ ALP) 

SEC2 = SEC * SFC 
SEC4 = SEC 2 * SFC 2 
SECfi = SEC A * SEC? 
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CSC=1.0/ SIN(ALP) 

CSC? = CSC * CSC 
CSC4 = CSC? * CSC2 
CSC6 = C SC 4 * CSC2 
CDT-WO/TN 

CF1=AZ*ALGG{ ( SEC+TN)/( SEC-TN) )*W*i-2.0) 

WA=W*AZ/S 
WA 2 = WA * WA 

WA4 = WA2 * WA 2 

WA 6 = Vi A4 * WA 2 

CF1=CF1*( 1.0-WA2 /12.0+WA4 /320.0- WA6 / 1 6 1 2 S *0 ) 

CF2=WA2 *W*AZ*TN*SEC/3.0 

CF?=CF2*{ 1.0-WA2 *ISEC2 +1 . 5 } /40. 0+WA4 +ISEC4 /3.0 + SFC? 

1 *F 5012 + F5D8 )/840.0) 

CF3=BZ*AL0G{ ( C SC+COT ) / ( C SC- COT) )*W* 1-2.0) 

WB~W*BZ /5 
WB2 = WB * WB 

WB4 = WR2 * WB2 

WB6 = WB4 * WB 2 

CF3=CF3*{ 1.0-WB2 /12.0+WB4 /320.0-WB6 /16128.0) 

CF4=WB 2 *W*BZ*CGT*CSC /3. 0 

CF4=CF4*C 1.0-WB2 *(CSC2 + 1- 5) /40. 0+WB4 *(CSC4 /3.0 + CSC2 

1*F5012 + F5D9I/840. 0) 

PTT=CF1 +CF2+CF3+CF4 
DF 1 = WA 2 *S*TN*2.0 

OF l=CFl*( 1,.0-WA? * { SEC 2 +2- 0) /36. 0+WA4 =M3.0*SEC4 +4.0*SFC? 

1+8.0 )/54C0„0-WAfi * [15* 0* SFC6 +18.0*SFC4 +24.0*SFC? +48*01/2116 
2800.0 

DF2=WB2 *S*C0T*2,Q 

nF2=DF2*t 1.0-WB2 *{CSC2 + 2. 0) /36. 0+WB4 *{3.0*C5C4 +4.0*CSC2 + 

18.0 } /5400. C-WB6 *{CSC6 +6.0*CSC4 /5. 0+8. 0*CSC 2 /5.0+F16D5 )/ 

2141120.0 
PBB=0*=1+CF2 
GO TO 2£1 
235 CONTINUF 

nn 250 J=1,LP 
ILK=J 

ILL ~NN + 2*{ J— 1 ) 

P B{ J 1=0. C 
PT( J ) = 0. C 

CALL GAUS S I X » WT» ILL) 

SANR =C ■ C 
SAN 1 = 0*0 

THESE LOOPS ARE FOR THF NUMERICAL INTEGRATION 

DO 255 II=1 T ILL 
R I =RR- X( II) * BZ 

DO 255 J J=1 » ILL 
T I = TT-DT* X ( J J ) / ?. 0 

R3=( (R*CT-RI* COSITI) ) **2+ l R* ST-R I * SI N ( TI )) **2 ) **0 . 5 
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ART=RO*WDS 
CART- CCS ( ART 1 
SAR T = S IN( ART ) 

S IN T=R I *W* CDS! FMW *TI) /RO 

SANR-SANR +SINT* SART*WTU I )*WT( J J ) 

SANI=SANI+$INT*CART*WT{T I)*WT(JJ)*(-1.0) 

255 CONTINUE 

PBl J )=SANR*B70T 

PT ( J)=SANI*B*0(K)*DT M. 0 

PT t J ) = SAN I + 8ZOT 

IF IJ .LT. 2) GO TO 250 

JL=J-1 

ABT= ABS( ABSIPBCJ) }- ABSIPB(JL)J) 

ARC- AB S { ABSIPTUn- ABSIPTUDJ) 

ARTC-ABT-CHK 

ABCC=ABC“CHK 

IF ( ABTC.LT. C-DO • AND • ABCC . L T. 0* DO ) GO TO 260 

257 IF(J-LP) 250,253,260 
253 CONTINUE 

WRITE (6,1701 TT,RR,ABT V ABC 

170 FORMAT ( 1H0, 1 INTEGRA TI ON FOR BLOCK CENTERED AT THET A= 1 ,1PF13*6, 

1* R = 1 , E13* 6 ,/,19X,'0JD NOT C ONVERGE * , / , 19 X , * RE AL DIFF =',E13.6 
2 , « IMAG OIFF = * ,F13.61 

250 CONTINUE 

260 CONTINUE 
PBB=PBt ILK) 

PTT^PTt ILK) 

261 CONTINUE 

PART= CMPLX I PBB t>PTT ) 

PO=PO+PART*PHtO(K) 

230 CONTINUE 
220 CONTINUE 
RETURN 
END 
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SUBROUTINE PRESHR) 


Z0M915 21 - 21 

THIS SUBROUTINE CALCULATES THF CONTR I RUT ION OF PHT1 
TO THE PRESSURE ON THE DUCT EXIT PLANE 
COMPLEX *8 CMPLX 

1 , PH 10(40} , PHI1 (40) ,PHI2 (40) ,P0»Pl,P2,PART 

DIMENS TON X(?0),WT(?0) ,0(40} t C(40) ,F(40} ,PB(40) ,PT(40) 
COMMON /DATA /l ,M,N,NN,LP , MW, A , B T W , S ,CHK T 0 ,C , E , P I 
COM MON/ PH I /PHI 0,P Hi 1, PH I 2, PC, PI ,P2 

COMMON /ANG/ FINK t BD2, WDS, FMW, PI2, RW, AW,BB, AA 
Pl= CMPLX(C.OFOtO.OEO) 

LOOP FOR THE TANGENTIAL VARIABLE 

CD 230 K=1 t N 
[F(K-1J221, 231,232 

231 ZZ=-C(K)*BD2 
GO TO 233 

232 KL-K-1 

ZZ=ZZ“(C(KU+C(K) )*BD2 

233 CONTINUE 

BCD? = B * C ( K ) / 2. DO 

LOOP FOR THE ANGULAR VARIABLE 

00 ?7C 1=1, L 
IF { 1-1 } 221,221,222 

2 21 TT=0.0 
GO TO 2 23 

222 IL = 1-1 

TT=TT+(E(IL )+E (I) )*PI 
??3 DT = E ( I ) * PI? 

BCDT4 = BCD2 * CT / 2.00 
DO 250 J=l,LP 
IL K = J 

ILL=NN + ?*{ J— 1 J 
CALL GAUSS( X, WT, ILL) 

PB( J ) = 0.C 
PT( J )=C.C 
SANR=O.C 
SAN I =0 ■ C 

THF SF LOOPS ARE FOR THE NUMERICAL INTEGRATION 
C 

no 255 11=1, ILL 
Z I = ZZ- XI III * BCD2 

DO 255 JJ=1,ILL 
Tr = TT-DT+XC JJJ/2.0 

RI = l (R-B* CDSITIJ )**?+( RR )*( SINtTI ))**2+ZI*ZT)**0.5 
ART=RI*fcCS 
C ART= CCS ( ART ) 
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SARI- S INt ART ) 

5INT=BW * COS! EMW *Tfi/RI 

SANR=SANR+S INT* SART*WT( II ) *WT (J J ) 

SANI-SANI+5 INT*CART*WT ( II )*WT(JJ)*{*“1«0) 

255 CONTINUE 

PB t J)-SANR*BCDT4 
PT t J)-SANI*BCDT4 

IF (J .LT. 2 ) GO TO 250 

JL=J-1 

ABT= AB SI ABStPBtJ) )- ABStPBtJUM 
ABO ABSt ABSlPTUIl- ABStPUJL))} 

ABTC=ABT-CHK 

ARCC=ABC-CHK 

IFU ABTC.LT.O. I ,AND. ( ABCC* LT. 0. ) ) GO TO 260 
257 IF(J-LP) 250,253,260 
253 CONTINUE 

WRITE (6,170) TT,ZZ ,ABT,ABC 

170 FORHATUHOi 'INTEGRATION FOR BLOCK CENTERED AT T HET A*» , 1 PE13*6, 
1* Z=* * E13* 6 ,/,19X, *DID NDT CONVERGF* ,/,19X,*REAL DIFF -»,E13 # 6 

2 , 1 IMAG DIFF = » ,E13.6) 

250 CONTINUE 
260 CONTINUE 

PBB=PB( ILK) 

PTT^PTI ILK) 

PART- CMPLX { PBB ,PTT } 

P l^P 1+P ART*PHl 1 (K } 

270 CONTINUE 
230 CONTINUE 
RETURN 
END 
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SUBROUTINE PRES2! R ) 

ZDM915 21 - 21 

THIS SUBROUTINE CALCULATES THE CONTRIBUTION OF PHI2 
TO THE PRESSURE ON THE DUCT EXIT PLANE 
COMPLEX *8 CMP LX 

COMPLEX *8 PH 10 { 403 t PH 1 1 (40) , PHI 2 (40) , PO , PI * P2 ,P ART 

DIMENSION X( 20 1 ,WT!20) » 0 t 40 1 ,C ( 40) , E (4-0 ) ,P B [40 ) T PT £ 40 ) , C? £40 } 
COMMON /PHI/PHI CtPHI 1, PH 1 2 , PO , PI * P2 
CnMMnN/rATA/L T M,N,NN,LP,MW,A,B,W,S,CHK,D,C,EtPI 
COMMON ZC2N2/C2,N? 

COMMON/ANG/ FINK, BD2, WDS, FMW, PI2 , BW, AW,BB, AA 
P2- CMPLXIQ.OEO,O.OEO) 

105 FORMAT! 1H ,2X,4D25. 14) 

LOOP FOR THE TANGFNTIAL VARIABLE 

00 230 K=1,N2 

TF t K- 1 J 52 1,531, 532 

531 ZZ=C2IK)*BD2 
GO TO 533 

532 KL =K- 1 

ZZ=ZZ+(C2(KL)+C2£K) )*BD2 

533 CONTINUE 

BCD2 = B *C 2£ K ) / 2.D0 

LOOP FOR THE ANGULAR VARIABLE 

‘30 27C 1 = 1, L 
lF(I-l) 2' 1 ,221 ,222 

? 2 i rr=o.o 
GO TO 223 

222 TL= I— 1 

TT =TT+ 1 F ( IL H-E £1) I* PI 

223 DT= E ( r )* P 12 

BDDT4 = B * D £ K J * OT / 4.D0 
WRITE{6,105) TT,ZZ 
DO 250 J= 1 , LP 
WRITE! 6,5) J 
5 FORMAT! IF ,10X, 15) 

ILK = J 

ILL =NN +2* ( J- 1 1 
P8 ( J ! = G »C 
PT( J J = 0*C 

CALL GALSS(X,WT t ILL) 

SANR =0 • C 
SAN 1 = 0 • C 

THFSE LOOPS ARE FOR THE NUMERICAL INTEGRATION 

DO 255 11=1, ILL 

ZI=ZZ- X£ II ) * BCH2 
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DO 255 JJ=l,rtL 
TI=TT-DT*X( JJ 1/2.0 

RI = UR-A* C0StTin**2+ AA *{ S I N{ TT ) ) **2+Z I *Z I ) **0*5 
ART^R I*fcCS 
CART= COStART) 

SART = SINtARTI 

S INT=AW * COS* FMW *TU/RI 

SANR=S ANR +S INT* SART*WT t 1 1) * WT U J ) 

SAN I ^SAN I+S INT*CART*WT III J * W7 t J J ) * t-1. 0 ) 

255 CONTINUE 

PB( J ) =S ANR+BDDT A 
PT (J l=SANI*BDDT4 

IF (J «LT« 2) GO TO 250 

JL = J— 1 

ABT = ABS{ ABSIPB ( J ) )~ ABStPStJUM 
ABC= ABSt ABSIPT(J) 1- ABStPTtJLHI 
ABTC = AB T-CHK 
ABCC=ABC-CHK 

TF[ ( ABTC.LT*0. 1 .AND . ( ABCC. LT. 0. )} GO TO 260 
IF { U-LP ) 250*253*260 
253 CONTINUE 

WRITE (6,170 TT , ZZ , A3T , ABC 

170 FORMAT t 1H0, 'INTEGRATION FOR BLOCK CENTERED AT THET A-» f lPE13.6, 
1* Z = ', E13.6 ,/,l<5X,'DID NOT CONVERGE *,/ ,1 9X, • REAL DIFF =',E13.6 

2 IMAG DIFF = *,E13.6> 

250 CONTINUE 
260 CONTINUE 

P8B=PB ( ILK) 

PTT=PT( ILK > 

WR ITEt 6 , 106 5 PBB ? P TT 

106 FORMAT! 1H ,10X,'P2 REAL = ',E25.7 ,« P? rMAG^',E25.7 /) 

P ART= CMPLX ( PBB ,PTT ) 

P2=P2*PART*PHI2(K) 

270 CONTINUE 
230 CONTINUE 
RETURN 
END 
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APPENDIX 7 

GAUSS SUBROUTINE LISTING 


^RECEDING RAGE BLANK NOT EJLMEQ’ 
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SUBROUTINE GAUSSt Xt WT t H) 
DIMFNS ION X I 1 2 ) »WT( 12 ) 

IF t M-A ) 10 » 1 1 r 1 2 

10 Xtl)=-0. 577350269189626 
Xt 2) = -Xt 1 1 

WTU) = 1.C 
WTt 2 ) = 1.C 
GO TO 20 

11 Xtl)=-0. 861136311594053 

X( 2) =-0.335981043584856 
X t 3 ) =- X t 2 1 

X ( 4 ) Xtll 
WT(l) =0.347854845137454 
WTt?) =0.6 52 1451 54862546 
WTt 3)=MT( 2) 

WTt 4)=WTt 1 ) 

GO TO 20 

12 IFtH-81 13,14, 15 

13 Xt l)=-0. 932469514203152 
Xt 2) =-0.6612093 86466265 
Xt 3) =-0.23861 9186083197 
X t 4 J =— Xt 3 ) 

X t 5 ) =- X t 2 ) 

X t 6 ) =- Xtl) 

WTt 1) =0.1 7 13244523791 70 
WTt 2 ) = C. 3607 615 73 0481 39 
WT(3 )=0. 467913934572691 
WT(4l=WTt 3» 

WTt 5)=WT(2 ) 

WTt 6)=WT( 1 ) 

GO TO 20 

14 Xtl )=-0. 960289856497536 
Xt 2) =-0.796666477413627 
Xt 3)=-0.5255324C99l6329 
X ( 4 ) =- 0 . 1834 34 642495650 

XI 5 )=- X t 4 ) 

X t 6 ) =- X t 3 ) 

X { 7 ) =— X ( 2 ) 

X l 8 ) =- XU) 

WTtl ) =0.10 122 85362903 76 
WTt 2 ) = 0. 22 238 10344533 74 
WTt 3 1=0.313706645877887 
WT(4)=0.362683783378362 
WTt 5 J=WT14) 

WT ( 6 )=WT ( 3 ) 

WTt 7)=WTt2) 

WTt 8 ) =WTC 1 ) 

GO TO 2C 

15 I FI M- 12 ) 16, 17, 18 

16 Xfl)=-0. 5739065285171 72 
X ( 2) =-0.865063366688985 



XI 3)=-0. 679409568299024 
XI 4> =-0.433395394 12 9247 
X( 5 )=— 0 . 148 8743 38 981631 

X ( 6 > -- X ( 5 ) 

X I 7 ) =- X{ 4) 

XI 8 >=- X I 3 ) 

X ( 9 } =— X ( 2 ) 

xi io i =- xm 

WT ( 1 1 = 0. C6 6671344308688 
WT[ 2 >=0.149451349 15 0581 
WTI 3 > = 0.219086362 51 59 82 
WT { 4 > = 0.269266719309996 
WT ( 5 >=0.295524224714753 
WTI 6 J=WT( 5} 

WTI 71=WT{4) 

WT{ 8> = WT( 3) 

WTI 9 I =WT I 2 ) 

WT( 1 0 1 = W T I 1> 

GO TO 20 

17 XI l)=-0. 98156063424671 9 
XI 2>=-0. 904117256370475 
XI 3) =-0.769902674194305 
XI 4 >=-0.58 731 79542 8661 7 
XI 5 >=-0.36783149891 8180 
XI 6 >=-0.1253334 08511469 

XI 7>=- X I 6 ) 

X(8)=- X ( 5 ) 

X ( 9 } =— X ( 4 ) 

XI 10 ) =- X ( 3 } 

XI1U=- X ( 2 ) 

xc i2>=- xm 

WTI 1 >=0.047175336386512 
WTI 2 > = 0.1069 393 2599 5318 
WTI3 > = 0.160078328543346 
WTI 4 >=0.203 1674 26723066 
WTI 5>=0. 233 49 2536 53 83 55 
WTI 6 >=0.24914704 5 81 3403 
WTI 7>=WTI6I 
WTI 8>=WT( 5> 

WTI 9>=WTI4) 

WTI 10 ) =WTI 3 ) 

WTI 11 > =WTI 2 > 

WTI 12I=WTQ> 

18 IF IM-lfc) 

19 XI 1 >=-0.986283808696812 
XI 2 >=-0.928434883663574 
XI 3 >=-0.827201315069765 
XI 4) =-0.687292904811685 
XI 5 i =-0.515248636358154 
XI 6 >=-0.319112368927890 


Xt 7 1 =- 0 * 108054948707344 
XI 81 =- Xt 71 
X t 9 ) X t 6 ) 

Xt 10 J = - X t 5 ) 

X { 1 1 1 — - X ( 4 ) 

Xt 12)=~ X( 3 1 

XI 13|=- X I 2 1 

xc 14 ) =~ xm 

WT I 1 1=0. 035 11 94 £033 1752 

WTI 21=0.080158087159760 

WTI 31=0.121518570687903 

WTI 4 1 = 0. 15720316715 8194 

WTI 5 1 = 0. 185538397477938 

WTI 6 }=0. 205 19 8463 72 1296 

WT{ 7} =0*2 15263 8534631 58 

WT ( 8 1 — W T I 71 

WT ( 9 1 = W TI 6 1 

WT t 10 1 =WT( 5 1 

WT [ 1 1 1 = WT I 4 1 

WT { 121 = WT I 3 1 

WT i 13 1 =WT( 2 1 

WT ( 14 ) = WT ( 1 I 

21 Xtl)=-0. 989400934991650 
X{ 2) =-0*944 575023073233 
XI 3) =-0 *86 56312 0? 38 7832 
XI 4) =-0.755404408355003 
XI 5) =-0.617876244402644 
XI 6) =-0*458016777657227 
XI 7) =-0.28 16035 50779? 59 
Xt 81 =-0.095012 5 C983 7637 
XI 9 ) =- X t 8 ) 

X ( 101- — X l 7 } 

XI 11 )=- X( 6 ) 

x t i2i=- xm 

Xt 131=- Xt 4 ) 

XI 14)=- X I 3 J 

Xt 15)=— X ( 2 ) 

xi i6i=- xm 

WTI 11=0.027152459411754 
WT( 21=0*062253523938648 
WT{ 3 1=0.095158511682493 
WT{ 41=0.124628971255534 
WT ( 5 1=0.149595988816577 
WTI 6 1=0. 169156519395003 
WT(7 1=0.182603415044924 
WT( 8 1=0.189450610455069 
WTt 9l=WTt 8} 

WT { 10)=WT( 7) 

WTI U1 = WT( 61 
WTt 12)=WTI 51 


GO TO 20 
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WT{ 13 } =WT{ 4 1 
WT( 14 ) =WT ( 3 ) 
WT [ 1 5 I— ViTC 2 ) 
WTt 16»=feT* l ) 
20 RETURN 
FND 



APPENDIX 8 
EXAMPLE PROBLEM 
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0016400000000C 23C 8 

l.C 343,0 343.0 .0001 .03 

1.21 

2.500000- 02 2.5Q000D-02 5.00000D-02 5. 000000-02 5.00000D-02 5.000000-02 

7.500000- 0? 7. 500Q0D-02 7.50000D-02 7.50000D-02 7.500000-02 7. 500000-0? 

1.000000-01 1.000000-01 1.00000D-01 1 . 00Q00D— 01 1.000000-01 l.OOOOOD— 0 l 

1.000000- 01 l.OOOOOD-Ol 1.000000-01 1,250000-01 1.250000-01 1.250000-01 
1.250000-01 1.250000-01 1.25000D-01 1.25000D-01 1. 250000-01 1. 250000-01 
U 250000-01 1.250000-01 1. 250000-01 1. 250000-01 1*250000-01 W 500000-01 
lc 500000-01 1 , 5000CD-0 1 1. 50000D-01 1.50000D-01 

l. 25000D-G1 1.250000-01 1. 000000-01 7.500000-02 5.00000D-02 2.500000-02 

2.500000- 02 5.00000D-02 5, OOOOOD-02 7.500000-0? 7.500000-02 7.500000-02 

5. 000000- 02 5. CC000D-02 2.500000-02 2.50000D-02 


1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 




-6.051150-02 4. 396280-02- 7.07061D-02 4. 862 960-02 -9. 074RRD-02 5.779730-02 
-1. 16465D-01 6.955110-02-1.4271 80-01 8. 1 5449D-02-1 * 62 7780-01 9.070810-0? 
-1.784120-01 5. 7 8 49 2D- 02- 2, 06 01 0D-01 1.104570-01-2.527120-01 1.317960-01 
-3.338090-01 1.688750-01-4.855940-01 2,383490-01-7.459170-01 3.576790-01 
-1.10192000 5. 209690-01- 1.360750 00 6. 38583D-01-1. 684620-01 8. 127090-02 
9.232320 CC-4. 292080 00 

-5. 246290 01 2.439520 01 3. 30565D 00-1.525480 00 9.59800D 00-4.440850 00 
6.309650 00-2.905790 00 3.517480 00-1.605960 00 1,942270 00-8.739000-01 
9.675640-01-4.221780-01 4. 5611 70-01- 1 . 866000-01 2.335750-01-8.533070-02 
1.247160-01- 3. 677370-02 6. 554990-02- 1 . 1 1 6980-02 3.098380-02 3.159380-03 
7.Q2919D-03 1.24469D-02- 8. 372440-03 1 . 76 862D-02 -1 . 633480-0? 1.96714D-0? 
-2. 006600-02 1.98723D-02- 2.134450-02 1 . 908 20D-0? -2 • 119 72D-02 1.774860-02 
-2.0 24560- C2 1.6 15P2D-02- 1, 8901 40-02 1.446720-02-1.740920-02 1.2B081D-02 
-1.570470-02 1. 101020-02- 1. 392930-02 9. 15983D-03-1 . 2361 10-02 7.506910-03 
-1. 101260-02 6. C4770D-03- 9. 86 77 3D- 03 4. 76274U-03-8. 895X80-03 3.630740-03 
-8.061170-03 2.641050-03- 7.33882D-03 1, 770590-03-6. 69152D-03 1.004170-03 
-6. 1 07 44D-03 3.3 11800-04- 5. 56 072D-03-2 . 6 87300-04-5 . 037840-03 - 7. 98 58 1 0-04 
-4. 531750-C3- 1.27 128D-03- 4. 03 1010-03- 1.701 69D-03 -3. 468 190-0 3-2. 1488 80-0 3 
-2. 8 2 5960-03-2.643200-03-2. 13 0030-03-3 • 2 306RD-03-1 . 27398D-03 -4. 26R 3 10-0 3 
3, 542770- 04- 8. 188340-03 
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TUBE END IMPEDANCE PROGRAM FOR UNF LANGE C ANNULAR DUCT 


A = 0.0 

C HK= 9.9999Q9E-05 

L ~ 30 M - 16 

1. 250000E-QI 
5.000000E-02 
5.C00Q00F-07 
5. OOOOOOF— 0 2 

OUTER DUCT WALL 

2* 500000E— 02 
5.000000E-02 
7 • 499999 E-0 2 
9. 999996 F-0? 
9.999996E-02 
9. 99999 6E-02 
1. 250000E-01 
1. 250000E-01 
1.250000E-01 
1. 500000E-01 

THETA BOX WIDTH 

3. 000000 E-0 2 
3* 34-482 8E-0 2 
3.344828E-02 
3. 344828E-02 
3.34482BE-02 
3. 34482 8 F-0 2 
3. 34482 8 E-02 
3. 3 44 82 8 E-02 


INPUT DATA 


B = 1. OOOCOOF 00 K = 3. 430000F 02 S = 3.430000E 02 
R HO = 1. 209F99E OC 


N a 4C NN = 

1 » 25COOOE- 01 
2. 500000E- 02 
7.499999F- 02 
5, CC00C0E-02 


9.999996E— 02 
2* 50OC0CE-C2 
7. 4999996-02 
2. 500C00E— 02 


0 INK = 0 

7. 499999F-0 2 
5.0OCQOOF-O2 
7. 4999 99F-0? 
2* 5000'OOE-02 


BOX LENGTHS 


2. 5000 C0E-02 
5. COOOOOE- 02 
7. 499999E- 02 
9. 899996E- 02 
9.S99996E-02 
1 • 250000E - 01 
1. 250000E- 01 
1 • 25000CE - 01 
1. 250000E- 01 
1. 500QCOE-01 


5. 000000E-02 
7.499999E-C2 
7. 489999F— 02 
9. 999996E-0? 
9. 999996E-02 

1.2 50000E-C1 

1 . 2 50000F-01 

1.2 50C00E-01 
1 • 2 50000F-01 
1. 500000F-01 


5. OOOOOOF-02 
7. 499999E— 02 
7. 499999 E-02 
9.999996F-02 
9. 999996F-02 
1. 250 OOOE-O 1 
1. 250 OOOE-O "i 
l. 25000QE— 0 1 
1. 5 00 000 E-0 l 
1 • 500000E-0 1 


ON DUCT FACF 


3.344828E- 02 
3# 344828E- 02 
3.344828E-02 
3. 344828E- 07 
3.344828E- 02 
3. 344828E- 02 
3 • 3448 28E- 0? 
3.344828E-02 


3. 344828E-02 
3.344828E-02 

3.344 828E-02 
3.344828E-02 
3.344828E-02 

3. 344 02 8E-02 

3. 344 82 8F— 02 


3. 344 828F-02 
3. 344 8 2 8 F-0 2 
3. 344828E-0? 
3. 344 8 28 F-0 2 
3. 344828E-02 
3. 3448 2 8 F— 0 ? 
3. 344828E-02 
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VELOCITY ON DUCT FACE 


l.CCOOOOE 

00 

o.c 

l.OOOOOOE 

CO 

0*0 

l.OOOOOOE 

00 

o.c 

l.OOOCOOE 

oo 

0.0 

1* 000000 E 

00 

o.c 

1 • OOOCOOE 

00 

0.0 

l.COOOOOF 

00 

o.c 

l.OOOOOOE 

00 

0.0 

1.C000006 

00 

o.c 

l.OOOCOOE 

00 

0. 0 

l.OOOOOOE 

00 

0.0 

l.OOOCOOE 

CO 

0.0 

1. OOOQOOF 

00 

O.G 

1. OCOCOOE 

00 

0.0 

l.OOOOOOE 

00 

o.c 

1. OOOOOOE 

00 

0.0 


SOURCE DISTRIBUTION ON DUCT FACE. 


-6.051150E-02 
-9.074879E-02 
-1. 427180E-Q l 
-1. 78 412 OF-Ol 
-2. 5 27 120E-0 1 
-4. 85594 0E-G1 
-1. 101919E 00 
“1.60462OF-O1 


4.3962 80F- 02 
5. 779730E- 02 
8.1 544 B8E - 02 
9.784919E-02 
L. 317959E-0I 
2.383490E- 01 
5.209690E-01 
8. 127087E-02 


-7.070607E-02 
- 1* L64650F-C1 
- 1 • 62 7780E-01 
-2 * 06C iOOE-Ol 
-3.338C90E-01 
-7. 459 170E-01 
-Is 360749F 00 
9. 232 32 OF 00 


4. 862960F-0? 
6.955105E-0? 
9.070808F-0? 
1* 104569E-01 
1.688750F-01 
3.576789F-01 
6. 385829E-01 
-4. 292080E 00 


SOURCE DISTRIBUTION ON OUTER DUCT WALL 


-5* 24628 9F 01 
9.5S8000E 00 
3.517480E 00 
9. 675640F-C 1 
2.335750F— 01 
6. 5 5498 5E-02 
7.029187E-03 
-1.63 348 0E-Q2 
-2. 134450 E-Q2 
-2.024560E-02 
-1. 740920F-02 
-1. 392930F— 02 
— 1. 101260 E-0 2 
-8. 895 177F-03 
-7* 33881 8E-03 
-6. 1 0 7 43 8E-Q3 
-5. 0 378 3 7E-03 
—4.0 3 101 0F-03 
-2. 825960E-03 
-1. 273980F-03 


2.439519E 01 
-4.440849E 00 
- 1* E05960E 00 
-4. 2217 80E- 01 
-8 • 533067E- 02 
-U 116980E-02 
1. 244690E- 02 
1.9671 40E- 02 
1 . 9Q8200E- 02 
1 . 61 5220E- 02 
1. 280810E-02 
9. 159R29E-03 
6. C47700E- 03 
3. 630740F-03 
1. 770590E- 03 
3. 311 799E- 04 
-7. 985B09E- 04 

- 1. 701690F- 03 

- 2. 643 200E- 03 
-4. 2683C7E- 03 


3. 305650E 00 
6.309649E 00 
1.942269F 00 
4.5611 70F—0L 
1 .247 16CE-01 
3. 09838CE— 02 

- fi* 372 437E-C3 
-2.0C6600F-C? 
-2. 119720E-C2 
-1. 890140E-02 
-1^57047CE— 02 
-1.2361L0F-02 
-9. 867 72 8F — 03 
-8. 061 167E-03 
-6.691 5L9F-03 
-5. 560718E-03 
-4*531 749F— 03 

- 3 .46 8 190E -03 
-2.1 3C03CF— 03 

3.542770F-C4 


-1.525479E 00 
-2 . 905789E 00 
-8* 739000E-0 1 
-1. 866000E-01 
-3. 6773 70E-0 2 

3. 159380F-m 

1# 7686 20E-02 
1.987? 30E-0 ? 
1.774860E-02 
1.446720F-02 
1. 1010-20E-02 
7. 506907E-03 
4.762739E-03 
2.641050E-03 
1 • 004 1 70E— 03 
-2.687299 E— 04 
-1* 271 2 80F-03 
-2* 148880F— 03 
-3. 230680E-03 
-B. 188337E-03 
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INTEGRATION FOR BLOCK CENTERED AT THETA- 1.993285E-01 R= 

DID NOT CONVERGE 

REAL DIFF = 2.205372F-C6 I HAG PI FF = 


INTEGRATION FOR BLOCK CFNTERED AT THE TA= 6. CP3843F 00 R = 


DID NCT CONVERGE 
REAL DIFF = 2.324581E-C6 
ANNULAR RING NO 1 RADIAL POSITION - 

LOCAL IMPEDANCE = 
PRESSURE = 
VELOCITY = 


IMAG DIFF = 
6.250000E-02 
8*533 75 9E 01 
8.533759E 01 
l.OOOOOOF 00 


INTEGRATION FOR BLOCK CENTERED AT THE TA = 6.083343E 00 R= 

DID NOT CONVERGE 




REAL 

DIFF = 

3.814697E- 

C6 

I MAG DI'FF 

= 

ANNULAR 

RING 

NO 

2 

RADIAL 

POSITION 

= 

1 . 875C00F- 

01 





LOCAL 

IMPEDANCE 

= 

8* 916052 E 

01 






PRESSURE 

= 

8 • 916052 F 

01 






VELOCI TY 

= 

1 . OCOOOOE 

00 

ANNULAR 

R ING 

NO 

3 

RAD IAL 

POSITION 

- 

3. OCOOOOE - 

01 





LOCAL 

IMPEDANCE 


9.621 284 E 

01 






PRESSURE 

= 

9.621284F 

01 






VELOCI TY 

- 

l.OOOOOOE 

00 

ANNULAR 

RING 

NO 

4 

RAD IAL 

POSITION 

- 

3. 874999F - 

■01 





LOCAL 

I mpedance 

= 

1 . 045224F 

02 






PRESSURE 

= 

1.045224F 

0? 






VELOCITY 

= 

l.OOOOOOE 

00 

ANNULAR 

R ING 

NO 

5 

RADIAL 

POSIT! ON 


4. 499999E- 

■01 





LOCAL 

I MPFDANCE 


1.123271E 

02 






PRESSURE 

= 

1 .123271F 

02 






VELOCI TY 

= 

1 . OOOOOOF 

00 


INTEGRATION FOR BLOCK CENTERED AT TH E TA - 0.0 


DID NCT CONVERGE 
RFAL DIFF = 6.6161I6E-C6 
ANNULAR RING NO 6 RADIAL POSITION - 

LOCAL 1 MP F DANCE = 
PRESSURE = 
VELOCITY = 


P = 


I MAG DI'FF = 
4.874998E-01 


1.179450E 
I.179450E 
1 .OOOOCOE 


02 

02 

00 


6.250000E-0? 

» 175781E-0? 

6.250000E-02 

5. 171 204F-02 

4.734263F 0? 
4# 734263F 02 
0.0 


1.875000E-01 

1.83L055F-04 

4. 804661E 0? 
4.B04661E 0? 
0.0 


4.938169E 0? 
4.936169E 02 
0 . 0 


5* 093 547 F 02 
5.093547E 02 
0.0 


5.243489F 02 
5.243489F 02 
0.0 


5. L24998E-C1 

1.983643F-04 

5.352490F 0? 
5.352490F 0? 
0.0 
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INTEGRATION FUR 

INTEGRATION FOR 
ANNULAR RING 

INTEGRATION FOR 
ANNULAR RING 

INTEGRATION FOR 

INTEGRATION FOR 
ANNULAR RING 

INTEGRATION FOR 
ANNULAR RING 


BLOCK CENTERED AT THE TA = 0.0 R= 

DID NOT CONVERGE 

REAL DIFF = 6.079674E-C6 IMAG DIFF = 


BLOCK CENTERED AT THE TA= 0.0 
DID NOT CONVERGE 
REAL DIFF = 5.722046E-C6 
NO 7 RADIAL POSITION = 
LOCAL IMPEDANCE = 
PRFSSURE = 
VELOCITY = 


R = 


IMAG DIFF = 
5.124998F-01 
1.223L649E 02 
1.271649F 02 
l.OOOOOOE 00 


BLOCK CENTERED AT THE TA^ 0.0 
DID NCT CONVERGE 
REAL DIFF = 1.621 246E-05 
NO 8 RADIAL POSITION = 
LOCAL IMPEDANCE = 
PRESSURE = 
VELOCITY = 


R = 


IMAG DIFF = 
5.499998E-01 
1.290967E 02 
1.290967E 02 
l.OOOOOOE 00 


BLOCK CENTERED AT THETA- 0. C R = 

DIO NOT CONVERGE 

REAL DIFF = 1.144409E-C5 IMAG DIFF = 


BLOCK CENTERED AT THETA= 6, CE3843E 00 R= 
OID NCT CONVERGE 

REAL DIFF = 1.811981E-C5 IMAG DIFF = 

NO 9 RADIAL POSITION = 5.999997E-01 

LOCAL IMPEDANCE = 1.399527E 02 

PRESSURE = 1.399527F 02 

VELOCITY = l.OOOOOOE 00 


BLOCK CENTERED AT THE TA= 6. C83B43F 00 R = 
DID NOT CONVERGE 

REAL DIFF = 1.525B79E-C5 IMAG DIFF = 

NO 10 RADIAL POSITION = 6.674997E-01 

LUCAL IMPEDANCE - 1.565026E 02 

PRESSURE = 1.565C26E 02 

VELOCITY = l.OOOOOOE 00 


4. 874998E-01 

2. 746 5 R2E-04 

5.499998E-01 

3. 356 934F-04 

5.4351 39 E 02 
5.435139E 02 
0.0 

6. 624997E— 0 1 

1.678467E-04 

5.571472E 0? 
5. 571 472 E 07 
0.0 

5.499998E— 01 

1. 220703F-04 

6.674997E-01 

2* 136230E-04 

5.787I5RE 0? 
5.7B7158F 0? 
0.0 

6.624997F-01 

1.983643F-04 

6.119973F 07 
6.119973E 02 
0.0 
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INTEGRATION FOR ELOCK CENTERED AT THETA = 0-0 K= 6 • 624997E-0 l 

DTD NCT CONVERGE 

REAL DIEF = 1.43C511E-C5 I HAG PIFF = 1.525R79E-04 

INTEGRATION FUR ELOCK C FNTERFD AT THE TA = 0- C R= 8.124996F-01 

DID NOT CONVERGE 

REAL DIFF = 2.574921E-05 IMAG DIFF = 4.272461E-04 

INTEGRATION FOR ELOCK CENTERED AT THETA = 6.C83843F 00 R= 7.374997E-01 

DID NCT CONVERGE 

REAL DIFF = 2.670288E-C5 IMAG D IFF « 7* 1362 30E-04 

INTEGRATION FOR BLOCK CENTERED AT THE TA = 6.C83843E 00 R= 8.124996E-01 

DID NGT CONVFRGE 

REAL DIFF - 2-384186E-05 IMAG DIFF = 2-28B8I8E-04 

ANNULAR RING NO 11 RADIAL POSITION = 7.374997F-01 

LOCAL IMPEDANCE = 1.830660E 02 6.667146F 0? 

PRESSURE = 1.830660E 02 6.667146E 02 

VELOCITY - l.OCOOOOF 00 0-0 


INTEGRATION FOR ELOCK CENTERED AT THE TA = 0. C R= 8.749996F-01 

DID NOT CONVFRGE 

REAL DIFF - 6.675720E-C6 IMAG DIFF = 7. 136230F-04 

INTEGRATION FOR BLOCK CENTERED AT THE TA - 6.C83843E 00 R= 7- 374997E-C1 

C ID NCT CONVERGE 

REAL DIFF = 1.049042E-05 I Ma G DIFF = 2-2888I8E-04 

ANNULAR RING NO 12 RADIAL POSITION = 8-124996E-01 

LOCAL IMPEDANCE = 2.218372E 02 7.467739E 0? 

PRESSURE = 2.218372E 02 7-467739F 0? 

VELOCITY = l.OCOOOOF 00 0.0 


INTEGRATION FOR BLOCK CENTERED AT THETA- 0.0 R= 8.12A996E-01 

DID NOT CONVERGE 

REAL DIFF = 1.335144E-C5 IMAG DIFF = 1.220703E-04 

INTEGRATION FOR ELOCK CENTERED AT TMETA= I.9932R5E-01 R= 0.124996E-O1 

DID NCT CONVERGE 

REAL DIFF - 1.621246E-05 IMAG U IFF = 1.983643F-04 

ANNULAR RING NO 13 RADIAL POSITION = 8.749996E-01 

LOCAL IMPEDANCE = 2.7363I1F 02 8. 5607R6F 02 
PRESSURE = 2. 73631 I E 0? 8.560786E 07 
VELOCITY = l.OOOdOOE 00 0.0 
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R = B.749096E-01 


INTEGRATION FOR BLOCK CENTERED AT THE TA-= 0.0 

DID NOT CONVERGE 
REAL DIFF = 1.335144E-C5 [MAG DIFF = 


INTEGRATION FOR 


ANNULAR RING 


BLOCK CENTERED AT THE TA = 0.0 
DID NCT CONVERGE 
REAL Cl FF - 9 • 536 743E- 06 
NO 14 RAD TAL POSITION * 
LOCAL IMPEDANCE = 
PRESSURE = 
VELOCITY = 


R = 


IMAG DUFF = 
9.249995E-01 
3.457480E 02 
3.457480E 02 
l.OOOOOOF 00 


INTEGRATION FOR BLOCK CENTERED AT THETA” 0.0 R= 

DIO NOT CONVERGE 

REAL DIFF = I.I44409E-C5 IMAG DIFF = 


INTEGRATION FOR ELOCK CENTERED AT THE TA = 0.0 


R = 


DID NCT CONVERGE 
REAL DIFF = 7.629395E-06 
ANNULAR RING NO 15 RADIAL POSITION = 

LOCAL IMPEDANCE = 
PRESSURE = 
VELOCtTY = 


IMAG DIFF = 
9.624995E-01 
4.480442E 02 
4.480442E 02 
I • OOOOOOE 00 


INTEGRATION FOR BLOCK CENTERED AT THETA* 0.0 R= 

DID NCT CONVERGE 

REAL DIFF - 6.675720E-C6 IMAG DIFF - 


INTEGRATION FOR BLOCK CENTERED AT THE TA= 0.0 l = 

DID NOT C ONVFRGE 

REAL DIFF = 4.768372E-C6 IMAG DIFF = 


INTEGRATION FOR 


ANNULAR RING 


BLOCK CENTERED AT THE TA = 0 
DID NCT CONVERGE 
REAL DIFF = 5.722046E-06 
NO 16 RADIAL POSITION = 
LOCAL IMPEDANCE * 
PRESSURE = 
VELOCITY = 


.0 

IMAG DIFF = 
9.874995E-01 
6.368638E 0? 
6.368638E 02 
1* OOOOOQF 00 


Z = 


RESISTANCE RATIO, TAU = 0.5I969945E 00 

REACTANCE RATIO, CHI = 0.19329624E 01 


1.831 055E-04 

9.624995E-0I 

2.899170F-04 

I. 0I0O26F 03 
1. 010026E 03 
0.0 

9.249995E-01 

2. 1 86 5 84E-02 

9.874995E-01 

2.427673E-02 

1.229231E 03 
I.22923IF 03 
0.0 

9.62 4995 F— 0 1 

2.4R5657E-02 

-1.250000E-02 

7. 09? 2 95 E-02 

—3* 75 0 000 E -02 

1.83I055E-04 

I* 635084E 03 
I. 6350R4E 03 
0.0 
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